import cv2
import numpy as np
import json


class MyPoint:
    def __init__(self, x , y , flag):
        self.point = (x, y)
        self.flag = flag





# 读取TIFF文件
tiff_file = '/Users/daxiang/Downloads/haoping.tif'
png_file = '/Users/daxiang/Downloads/out.png'
image = cv2.imread(png_file, cv2.IMREAD_GRAYSCALE)

has_black_border = False
# 设置缩放比例
scale_percent = 1  # 缩放百分比，这里设置为10%

# 添加50个像素的黑边: todo 添加黑边，是右边，去掉黑边是左边
if has_black_border:
    border_color = (0, 0, 0)  # 黑色
    border_width = 50
    image = cv2.copyMakeBorder(image, border_width, border_width, border_width, border_width, cv2.BORDER_CONSTANT, value=border_color)


for i in range(2):
    image = cv2.GaussianBlur(image, (3, 3), 0)

# 阈值分割，将图像转换为二值图像
_, binary_image = cv2.threshold(image, 0, 255, cv2.THRESH_BINARY)

#
kernel = np.ones((3, 3), np.uint8)
# 膨胀操作去除内部杂点噪音  eroded
for i in range(3):
    binary_image = cv2.dilate(binary_image , kernel, iterations=1)
    # 腐蚀操作还原轮廓形状
    binary_image = cv2.erode(binary_image, kernel, iterations=1)


# cv2.imshow('Image', binary_image)
# cv2.waitKey(0)
# cv2.destroyAllWindows()


# -------------  dst = cv2.Canny(binary_image,10,50)
dst = cv2.Laplacian(binary_image,-1,ksize=3)

# # # 检测轮廓
contours, _ = cv2.findContours(dst, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

# # 轮廓近似
epsilon = 0.01 * cv2.arcLength(contours[0], True)
approx = cv2.approxPolyDP(contours[0], epsilon, True)

# 获取地理转换信息（根据您的TIFF文件的实际情况进行调整）
geotransform = (111.2279051057111, 0.000011751600569, 0, 32.85772767224697, 0, -0.000009912427693)

contour_image = np.zeros_like(image)
cv2.drawContours(contour_image, contours, -1, 255, 1)

# 显示图片
cv2.imshow('Image', contour_image)
cv2.waitKey(0)
cv2.destroyAllWindows()



# # 开始提取边界坐标
# # 提取轮廓边界线上的所有点的坐标
points = np.argwhere(contour_image == 255 )
center = np.mean(points, axis=0)
# leftmost_point = max(points, key=lambda p0: p0[0])
coordinates_res = sorted(points, key=lambda p: np.arctan2(p[1]-center[1] , p[0]-center[0]-30))

print(f"------len(points): ", len(points))
print(f"------len(coordinates_res): ", len(coordinates_res))
print(f"------coordinates_res: ", coordinates_res)
#
new_coordinates_res = list()
for index in range(0,len(coordinates_res), 2):
    new_coordinates_res.append(coordinates_res[index])

print(f"------len(new_coordinates_res): ", len(new_coordinates_res))
new_coordinates_res = np.array(new_coordinates_res, np.int32)
new_coordinates_res = new_coordinates_res.reshape((-1, 1, 2))

new_contour_image = np.zeros((269, 229, 3), dtype=np.uint8)

cv2.polylines(new_contour_image, [new_coordinates_res], isClosed=False, color=(255, 255, 255), thickness=1)
cv2.imshow('Image', new_contour_image)
cv2.waitKey(0)
cv2.destroyAllWindows()








latitudes = list()
longitudes = list()
for index in range(len(coordinates_res)):
    point = coordinates_res[index]
    # 转换为经纬度坐标： todo 带黑边，坐标需要减44，
    if has_black_border:
        latitudes.append(geotransform[3] + (point[0]-44) * geotransform[5] * (100 / scale_percent))
        longitudes.append(geotransform[0] + (point[1]-44) * geotransform[1] * (100 / scale_percent))
    else:
        latitudes.append(geotransform[3] + (point[0]) * geotransform[5] * (100 / scale_percent))
        longitudes.append(geotransform[0] + (point[1]) * geotransform[1] * (100 / scale_percent))


# 创建GeoJSON对象
geojson = {
    "type": "Feature",
    "geometry": {
        "type": "LineString",
        "coordinates": []
    }
}

# 添加经纬度坐标到GeoJSON对象
for latitude, longitude in zip(latitudes, longitudes):
    geojson["geometry"]["coordinates"].append([longitude, latitude])

# 将GeoJSON对象转换为字符串
geojson_str = json.dumps(geojson)

output_file = '/Users/daxiang/Downloads/haoping.json'
with open(output_file, 'w') as f:
    f.write(geojson_str)



